Enhanced optimization-based method for the generation of patient-specific models of Purkinje networks

Cardiac Purkinje networks are a fundamental part of the conduction system and are known to initiate a variety of cardiac arrhythmias. However, patient-specific modeling of Purkinje networks remains a challenge due to their high morphological complexity. This work presents a novel method based on optimization principles for the generation of Purkinje networks that combines geometric and activation accuracy in branch size, bifurcation angles, and Purkinje-ventricular-junction activation times. Three biventricular meshes with increasing levels of complexity are used to evaluate the performance of our approach. Purkinje-tissue coupled monodomain simulations are executed to evaluate the generated networks in a realistic scenario using the most recent Purkinje/ventricular human cellular models and physiological values for the Purkinje-ventricular-junction characteristic delay. The results demonstrate that the new method can generate patient-specific Purkinje networks with controlled morphological metrics and specified local activation times at the Purkinje-ventricular junctions.

. Illustration showing the different geometrical concepts related to the method alongside the data structure used to store the nodes and segments of the PN. In panel (A), the concept of a segment and branch is highlighted in blue and red, respectively. In addition, the initial root branch of the PN is colored in black, while an example of a terminal branch is presented in green. In panel (B), the data structures used to store a PN are illustrated. There are two arrays, one of nodes and another of segments which keep track of the current state of the PN. Each segment has access to its parent and left/right off-springs, which are all pointers to the segment structure. Moreover, each segment has two pointers to the two nodes that define the segment. Using these pointers is possible to retrieve the proximal and distal coordinates of any segment.

A.2.1 Biventricular surface
The simplified biventricular mesh was constructed by taking the differences between three ellipsoids, one for the Left ventricle (LV) and Right ventricle (RV) endocardium and one for the epicardium. The following equation gives an ellipsoid that is not centered on the domain origin: where c x , c y and c z are the center coordinates. The LV surface is built by taking the difference between two ellipsoids. The LV epicardium given by (3) is subtracted from the LV endocardium given by (2). Similarly, the RV surface is built by taking the difference between three ellipsoids. The RV epicardium given by (5) is subtracted from the RV endocardium (4) and the LV surface previously calculated. Next, the biventricular surface is assembled by joining the LV and RV surfaces and cutting the base plane defined by x = 0. Finally, the 1/23 mesh is translated to the origin and scaled to micrometers.

A.2.2 Purkinje network
To build the reference Purkinje network (PN) for the simplified mesh, we use the fractal algorithm from 1 with the parameters given by Table S1. Regarding the parameter description from the fractal algorithm by Costabal et al. 1 , the Root and Second node coordinates are utilized to give the initial growth direction of the first segment. Fascicle angle and length are related to the generation of the first segment by providing information if there will be a bifurcation after the first segment using a specified angle. If the Fascicle angle is zero, a linear fascicle segment grows with a designated length from the distal position of the first segment. The Repulsive parameter regulates the branch curvature: the larger the repulsion parameter, the more the branches repel each other. The branch angle controls the mean branch bifurcation angles. The number of iterations configures the number of growth iterations. After the networks were generated, the structures were translated to the origin and scaled to micrometers.  10 10 Table S1 . Parameters used for the generation of the reference PN of the Simplified mesh.

A.3 Purkinje-Ventricular-Junction delay tuning
To model the Purkinje-Ventricular-Junctions (PVJs) sites, an additional current, I tiss PV J , was included on the tissue cells that are coupled to the terminal cells from the PN as described in the following equation: where V purk is the transmembrane potential of the terminal Purkinje cells, V tiss i is the transmembrane potential of the ventricular cells i attached to the Purkinje cell, R PV J is a fixed-resistance and N PV J is the maximum number of ventricular cells that are inside the PVJ site. This additional current is included on the right hand side of the associate linear system of the ventricular domain only for the ventricular cells related to the coupling.
Similarly, for the Purkinje cells a current, I purk PV J , coming from the ventricular cells that are coupled to the terminal Purkinje cell, is included on the right hand side of the associate linear system of the Purkinje domain only for the terminal Purkinje cell related to the coupling and is given by the following equation: In Figure S2, we illustrate how the Purkinje coupling was implemented and in Figure S3 how the experiment was done. Figure S2 . Illustration of the Purkinje coupling model. The terminal Purkinje cells are linked to the nearest N PV J ventricular cells by a fixed resistance R PV J and two additional currents I purk PV J and I tiss PV J are computed following equations (6) and (7) for only the Purkinje and ventricular cells that are related to the coupling. Figure S3 . Experiment to calibrate the Purkinje-Ventricular-Junction delay in the anterograde direction. A single 10cm Purkinje fiber stimulates a face of a tissue block (1cm 3 ). The PVJ delay in the anterograde direction is measured by the difference between the LAT of the ventricular cells and Purkinje cell related to the coupling. The Purkinje cells are modeled using the Trovato human Purkinje model 2 and the ventricular cells using the recent version of the Tomek human ventricular model 3 . The action potential traces from both the Purkinje and one of the ventricular cells of the coupling are depicted alongside the PVJ delay when R PV J = 900kΩ and N PV J = 60 in the right panel.

2/23
Applying different combinations of R PV J , and N PV J we measured the PVJ delay as Table S2 illustrates. Based on the results from Table S2, the Purkinje coupling parameters were set to R PV J = 700kΩ and N PV J = 60, generating a PVJ delay of approximately 3 ms in the anterograde direction, which is within the acceptable range 4 .

A.4 Geometric results
In the following section the geometrical results from the 10 best PNs regarding the three biventricular meshes are presented. To evaluate the PNs geometrically the branch size (millimeters), bifurcation angle ( o ) and the total number of branches and angles in the trees are calculated considering the mean and standard deviation of its values.

A.4.1 Simplified mesh
In Table S3 the geometrical results for the Simplified mesh are illustrated. As can be seen in this table, the values in terms of branch size demonstrate that the RV trees have a longer size when compared to the LV ones. When the average and standard deviation values are compared to the reference PN, we can noticed that the values of branch size for the 10 best networks can vary depending on the endocardium region where trees were generated. For the bifurcation angle measurement in this mesh, we could verify that the average values from both biventricular regions were close to the reference values and that the angles for the RV trees have more acute shape. Finally, for the total number of branches and bifurcations in the trees the values were exactly the same as the reference ones, which was already expected since all the generated PNs have the same number of terminals as the reference. This can also be verified by the fact that all the active PVJs are connected to the minimum network.  Table S4 . Results for the geometric features from the best 10 Purkinje networks of the Canine mesh. In the case of the Canine mesh, the Reference network was provided from the works 5, 6 .

A.4.3 Patient-specific mesh
In Table S5 the results for the patient-specific could not be compared to a reference PN since only the location and LAT of the active PVJs sites were available. However, important topological observations can be made from the generated PNs. In terms of branch size, the LV and RV trees for this particular mesh present the highest values. This can be justified due to the fact that the Patient-specific mesh have the largest endocardium volume when it is compared to the Simplified and Canine meshes leading to prolonged branches in either biventricular regions.  Table S5 . Results for the geometric features from the best 10 Purkinje networks of the Patient-specific mesh.

A.5 Activation results
In the next section the activation results from the 10 best PNs of the three biventricular meshes are depicted. To evaluate the PNs electrically the minimum and maximum LAT (milliseconds), maximum LAT error (milliseconds), RMSE (milliseconds), RRMSE (%), ε < 2ms (%) and ε < 5ms (%) are evaluated at the active PVJ sites in the trees.

A.5.1 Simplified mesh
In Table S6 the activation results for the Simplified mesh are illustrated. As can be seen in the table, the generated PNs presented a minimum and maximum LAT close to the reference values in both biventricular regions. When the maximum LAT error is analyzed the trees presented a good accuracy with errors below 3ms. This observation is also sustained by RMSE error results below 1ms and almost the majority of the active PVJs being connected within a LAT error tolerance below 2ms.

A.5.2 Canine mesh
In Table S7 the activation results for the Canine mesh are highlighted. Based on these results, for the minimum and maximum LAT, the generated PNs demonstrate a good accuracy with values close to the reference ones in both biventricular regions. In terms of maximum LAT error, we noticed that the LV trees have the highest values when compared to the RV trees which can be justified due to the fact that are more active PVJs in the LV region of this mesh than the RV. Analyzing the RMSE, RRMSE, ε < 2ms and ε < 5ms metrics the PNs from both regions presented good accuracy with errors equivalent to the results from the Simplified mesh.

A.5.3 Patient-specific mesh
In Table S8 the activation results for the Patient-specific mesh are presented. Similarly to the previous meshes, the generated PNs performed well in terms of minimum and maximum LAT leading to values close to the reference values. Only in the RV region of the Patient-specific mesh we found that a particular PVJ could not be connected with a LAT error below 5ms. When the other electrical measurements are analyzed, we noticed that the LV trees presented a better approximation to the reference LAT at the active PVJ sites. In addition, an interesting observation that can be made for the RV trees is that the 10 best networks for this region presented a very similar LAT even with a very distinctly topology.

A.6.1 Simplified mesh
In Figure S4 the results of the monodomain simulations for the Simplified mesh are depicted. Firstly, in Figure S4A, it is illustrated the baseline LAT calculated when the reference PN is coupled to the biventricular tissue. Secondly, in Figures S4B  and S4C, it is shown the LAT obtained when the best and worst PN generated by the method are coupled to the biventricular tissue, respectively. Finally, in Figure S4D and S4E, it is presented the absolute difference between the baseline and approximate LAT given by the PNs, where the results for best PN is on panel D and the worst PN on panel E. From the results of Figure S4, we noticed that both the best and worst PNs of our comparison set show a feasible LAT map for the Simplified mesh. As can be seen in panels D and E of Figure S4, the best and worst PNs present a maximum absolute LAT difference around 3ms, which is acceptable in terms of LAT matching. Figure S4 . Results for the LAT maps when a coupled monodomain simulation is executed over the Simplified mesh. First, in the top panel (A), the active PVJs are activated following the LAT provided by the reference PN. Next, in the middle panels (B) and (C), the LAT maps of the best and worst PNs generated by our method are depicted, respectively. Finally, in the bottom panels (D) and (E), the absolute difference between the LAT maps from the best and worst PNs is calculated, respectively.

A.6.2 Canine mesh
In Figure S5 the results of the monodomain simulations for the Canine mesh are illustrated. Initially, in Figure S5A, it is depicted the baseline LAT calculated when the reference PN is coupled to the biventricular tissue. Next, in Figures S5B and  S5C, it is shown the LAT obtained when the best and worst PNs generated by our method are coupled to the biventricular tissue, respectively. Lastly, in Figures S5D and S5E, it is presented the absolute difference between the baseline and approximate LAT given by the PNs, where the results for best PN is on panel D and the worst PN on panel E.
From the LAT and absolute difference maps in Figure S5, we can conclude that the best PN was able to activate the ventricular tissue similar to the reference LAT in almost all the biventricular tissue, except in a specific region of the RV. On the other hand, the worst PN could not activate several PVJs at the correct LAT even with a good approximation when the cable equation was used. This behaviour could be related to the PVJ coupling parameters used in the monodomain simulation. Mainly, from these results, we can conclude that certain PVJs must activated with a LAT very close to each other to stimulate these regions. In addition, as can be seen in these results, the PVJ delay plays an important role in ventricular activation and depending on certain conditions might be responsible for different LAT. Figure S5 . Results for the LAT maps when a coupled monodomain simulation is executed over the Canine mesh. First, in the top panel (A), the active PVJs are activated following the LAT provided by the reference PN. Next, in the middle panels (B) and (C), the LAT maps of the best and worst PNs generated by our method are depicted, respectively. Finally, in the bottom panels (D) and (E), the absolute difference between the LAT maps from the best and worst PNs is calculated, respectively.

A.6.3 Patient-specific mesh -Isolating the effects of the PNs
To further investigate the differences between the two simulated LAT maps generated by the best and worst PNs in the patient-specific mesh, which is depicted in Figure 6 from the main manuscript, we consider in Figure S6 the difference between the two respective LAT maps to verify if the cause of the errors are associated to the ventricular propagation model since the effects of the PN will be isolated.
As can be seen in Figure S6 we show the LAT maps generated when using the best and worst Purkinje networks for the patient-specific mesh together with the absolute LAT difference calculated between the two maps. As can be see in this figure and also in Figure 6 from the main manuscript, some effects of the Purkinje are isolated, like the region around a PVJ that could not be activated properly in the RV by both networks. However, from Figure S6 we can noticed that the activation of the LV is different between the two networks around the apex region with an absolute difference around 3ms and 4ms, which confirms that the errors are not only caused by the ventricular propagation model but also by the actual difference in the LAT at the PVJ sites of the two PNs which is sustain by the fact that the topology of the two PNs are different. Figure S6 . Results for the LAT maps for the best (panel A), and worst (panel B), PNs generated by our method considering the patient-specific mesh alongside the absolute difference LAT map (panel C) calculated by the difference between the two LAT maps. The topology of the PNs is illustrated in black color.

A.7 Method description
The detailed main program of the method depicted in Algorithm 1, starts by receiving as an input: a set of points S with locations for the terminal branches and an additional set of active PVJ points, an intermediate points cost function CF i , an active PVJ points cost function CF a , maximum number of feasible segments to be considered in the intermediate point cost function N i , maximum number of feasible segments to be considered in the active cost function N a , connection rate constant L rate , LAT error tolerance for an active PVJ connection L error , characteristic length of the domain l d , proximal location of the root branch x prox and an optional initial PN to be considered as starting root. The main parameters of the method are detailed in Table S9.

Algorithm 1: Main program.
Data: S, x prox , [initial PN]. Input parameters: CF i , CF a , N i , N a , l d , L rate , L error Result: Purkinje network generated within the set of points S. 1 S i ← PreProcessing(S, x prox ); # Algorithm 2 2 S a ← ExtractActivePVJ(S); # Algorithm 3 3 k term ← RootPlacement(S i , l d , x prox , [initial PN]); # Algorithm 4 4 while (not pass one time through S i ) do 5 x term ← GenerateTerminal(S i , N i , CF i , l d , k term ); # Algorithm 8 6 Include x term in the tree; 7 k term ← k term + 1; 8 if (k term % L rate == 0) then 9 AttemptPVJConnection(S a , N a , CF a , L error , l d , k term ); # Algorithm 10 10 end 11 end 12 PostProcessing(S a , N a , CF a , L error , l d , k term ); # Algorithm 13 13 Compute metrics and save network topology to a file;  Table S9 . Main parameters of the method.

A.7.1 Pre-processing steps
The PreProcessing subroutine described in Algorithm 2 and illustrated in Figure  The ExtractActivePVJ subroutine presented in Algorithm 3 and illustrated in Figure S7B, manages the set S by extracting the points that are considered active PVJs. In particular, each active PVJ must have an associated LAT value as a reference for comparison. With this value, we sort the active PVJ points using the QuickSort algorithm 7 .

A.7.2 Root placement
The RootPlacement subroutine detailed in Algorithm 4, is responsible for generating or loading the initial root of the PN. The method starts by constructing the root branch, in which, given the proximal location x prox , a distal location x dist is randomly selected from the intermediate point set S i and must attend a distance criterion. After selecting a feasible position, a geodesic pathway which connects x prox to x dist is constructed. Alternatively, the user can also provide an initial root structure described 15/23 in a PN topology format file that can be used to initialize the root of the method. Geodesic pathway: Build a geodesic path from x prox to x dist 20 k term ← 1; 21 end 22 return k term ; The subroutine SortPoint described in Algorithm 5 is responsible for selecting a point from the intermediate point set S i . The points from the intermediate point set are selected sequentially, and to allow variability in the PNs, this selection is made using different gap values. This gap value is defined as the difference between the previous and current index selected in S i . This value is calculated by taking the modulo between a random integer number, computed using the default rand function from the C++ standard library, and a fixed offset parameter, which was defined as 4 in this work. To avoid an index out of bounds error, we always apply a modulo operation between the next index and the total number of intermediate points in set S i . An important aspect of this function is that when the next_id variable returns to the start of the S i array, the stop criterion flag of the method's main loop is activated.

Algorithm 5: SortPoint subroutine.
Data: S i , global prev_id Result: Index of the next point to be selected in the set S i .

A.7.3 Distance criterion
The subroutine DistanceCriterion depicted in Algorithm 6, controls if all segments in set P s attend the distance criterion for a given x term position.

Algorithm 6: DistanceCriterion subroutine.
Data: x term , d thresh , P s Result: Return True if all segments in set P s attend the criterion or False otherwise.  Figure S8, generates a new intermediate terminal branch to the PN. The generation proceeds by sorting a point x term from set S i using the sort point function presented in Algorithm 5. The prospective location x term is accepted if x term satisfies the distance criterion, which is defined by Algorithm 6, and if there is otherwise we select a new point from S i . In case the above procedure fails 10 times, i.e., d crit < d thresh , the threshold distance d thresh is decreased by a factor 0.95. This is repeated until the acceptance of x term .
Next, when x term has been accepted as the distal end of a new terminal branch, the nearest N i feasible segments to generate a new branch are computed and stored in set F s using subroutine FillFeasibleSegmentsIntermediate which is detailed in Algorithm 9. After this procedure, x term is temporarily connected to the nearest N i feasible segments in F s in the model by constructing a geodesic pathway that links the distal point of the candidate segment, x dist , and x term . For each candidate connection, the value of the cost function given by CF i is obtained using Algorithm 7. Finally, the candidate connection with the lowest value of CF p and whose connection does not generate a collision to other segments in the tree is adopted and made 17/23 permanent. In case all feasible segments cause collisions, we sort a different x term from S i and recheck the distance criterion.  The AttemptGeneratePVJ subroutine described in Algorithm 11 and illustrated by Figure S9, respectively, attempt to connect an active PVJ, x PV J , by generating a new terminal branch. The branch is only made permanent if the absolute LAT error between the best candidate branch and the reference LAT value of the target active PVJ is less or equal to L error . Otherwise, the PVJ returns to set S a and the entire branch is pruned.
Data: x PV J , N a , CF a , L error , l d , k term Result: Try to connect the x PV J within the LAT error tolerance with a new branch. The FillFeasibleSegmentsActive subroutine as presented in Algorithm 12 is responsible for populating the F s set for the active cost function with the nearest N a segments to the x PV J position. This is done by passing through all the segments in the current PN and computing an approximation of the LAT error to the x PV J reference value. This approximation is calculated by

19/23
the sum of the current LAT of a given segment s i and the LAT given by the line that links the middle position x M of segment s i to x PV J using the cable equation.
Data: x PV J , N a , F s , P s Result: Fill the feasible segments set F s to connect x term .
Append (i, error) to F s ; 6 end 7 Apply QuickSort algorithm on F s using error as comparison parameter; 8  Best evaluation x PVJ x dist x PVJ (A) (B) (C) Figure S9 . Illustration of the AttemptGeneratePVJ subroutine described in Algorithm 11. In panel (A), a prospective location for a new PVJ branch, x PV J , is selected after attending the distance criterion. Next, in panel (B), N a feasible segments are evaluated by the cost function CF a and ranked in a table. After the evaluation step, segment 3 is unfeasible for the connection since it generates a collision in the PN and segment 4 returns the minimum value for the cost function and is within the LAT error tolerance L error = 2ms, for that reason is considered the best evaluation. Finally, in panel (C), a new terminal PVJ branch is constructed by linking the distal position of segment 4, x dist , to the location x PV J .

A.7.7 Post-processing steps
In the case there are remaining active PVJs to be connected after the main loop, we apply a post-processing step described and illustrated in Algorithm 13 and Figure S10, respectively. The first step in this procedure is to prune all intermediate branches of the tree. A branch is considered to be intermediate if it is not directly part of a pathway that links an active PVJ to the root. Next, the LAT error tolerance constraint is dropped by setting L error = ∞, and we attempt to connect all unconnected PVJs in S a using the pruned tree. The first N a feasible segments sorted by the LAT error are evaluated using cost function CF a and the segment with the best evaluation is connected using a geodesic pathway to x PV J . However, after this procedure some PVJs still could not be connected due to the distance criterion. This scenario could happen if x PV J is already too close to the current tree. In this particular case, we drop the distance criterion for these PVJs and attempt to connect the point with a feasible segment which returns the minimum LAT error using a geodesic pathway. If there are no geodesic pathway possible for x PV J , we consider the 5 closest segments by distance to x PV J and force its connection to 20/23 the one that returns the minimum LAT error using a straight line, regardless of any restriction. After this step, all the active PVJs that the user-specified are connected and the geometric and activation metrics are computed in this topology, which is referred to as minimum network. if (x PV J is not connected) then 6 AttemptGeneratePVJ(x PV J , N a , CF a , ∞, l d , k term ); # Algorithm 11 7 if (x PV J is not connected) then 8 Remove the distance criterion for x PV J ;

A.8 Finite Volume Method for solving the monodomain model in Purkinje-coupled simulations
The monodomain model is represented by a reaction-diffusion equation given by: The Finite Volume Method (FVM) can be used to solve the monodomain model. A common technique that solvers of cardiac electrophysiology perform when solving the monodomain model is to use an operator splitting to divide the equation in two distinct problems: a parabolic PDE associated to the diffusion part and a non-linear system of ODEs related to the reaction part of the model [8][9][10] .
In this work we solve the monodomain model for two different domains: the Purkinje and Ventricular. Within this context, whenever we are considering the transmembrane potential for the Purkinje domain we denote, V PK , and for the Ventricular domain, V V T . In summary, to advance the simulation by one timestep we need to solve four problems, two linear system associated to the diffusion part of the Purkinje and Ventricular domains, and two non-linear system of ODEs related to the reaction part of the Purkinje and Ventricular cellular models.
Additionally, it is important to remember that the Purkinje and Ventricular domains are coupled at the PVJ sites. For that reason we need to consider two extra currents: I purk PV J and I tiss PV J . These two currents describe the fluxes coming from the Ventricular to the Purkinje domains, and from the Purkinje to the Ventricular ones, respectively, are given by equations (6) and (7). In the following equations we describe the reaction and diffusion parts of the monodomain model for both the Purkinje and Ventricular domains after applying a time discretization using an implicit Euler scheme.
Firstly, we solve the reaction part of the Purkinje cellular model associated to equations (9) and (10), then we update the right-hand side of the Purkinje linear system which is related to the diffusion part of the model given by (11). Then, we compute the PVJ coupling current from the tissue to the Purkinje by adding the I purk PV J current coming from the coupled ventricular cells to the terminal Purkinje cells using equation (7), updating their transmembrane potential. During this step the current I purk PV J is added together with the fluxes from the faces of a control volume when the term ∇.(σ ∇V n+1/2 PK ) in equation (11) is discretized using the FVM. Next, we solve the Purkinje diffusion linear system given by (11). C m V n+1/4 PK −V n PK ∆t = −I ion PK (V n PK , η n PK ), ∂ η PK ∂t = f PK (V n PK , η n PK ), where V PK is the transmembrane potential for the Purkinje domain, I ion PK is the total ionic current for the Purkinje cellular model that may also depend on Purkinje gating variables η PK , I purk PV J is the PVJ coupling current coming from the coupled ventricular cells to the terminal Purkinje cells, λ is a binary variable that is equal to zero when a ventricular cell is not coupled to a Purkinje cell and one otherwise, σ PK is the monodomain Purkinje conductivity, β is the surface-volume ratio and C m is the membrane capacitance.
Secondly, we solve the Ventricular reaction part given by equations (12) and (13). Next, we update the right-hand side of the Ventricular linear system. Next, we compute the PVJ coupling current from the Purkinje to the tissue by adding the I tiss PV J current coming from the terminal Purkinje cells to the coupled ventricular cells using equation (6), updating their transmembrane potentials. Similarly, the current I purk PV J is added together with the fluxes from the faces of a control volume when the term ∇.(σ ∇V n+1 V T ) in equation (14) is discretized using the FVM. Finally, we solve the Ventricular diffusion linear system associated to equation (14) advancing the simulation by one timestep.

22/23
where V V T is the transmembrane potential for the Ventricular domain, I ion V T is the total ionic current for the Ventricular cellular model that may also depend on Ventricular gating variables η V T , I tiss PV J is the PVJ coupling current coming from the terminal Purkinje cells to the coupled ventricular cells, α is a binary variable that is equal to zero when a Purkinje cell is not coupled to a ventricular cell and one otherwise and σ V T is the monodomain Ventricular conductivity tensor.